home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_9 / a9_h.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  2.0 KB  |  77 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 9.S (Runge-Kutta Method for a higher order).
  9. % Section    9.7,    Runge-Kutta Methods, Page 475
  10. echo on; clc; format long; hold off; clear
  11.  
  12. % This program implements the Runge-Kutta method
  13. % for solving a second-order initial value problem:
  14. %     x'' = f(t,x,x')     with     x(a) = xa  
  15. %                                 x'(a) = ya
  16.  
  17. % Which is equivalent to solving the system:
  18. %      x' = y             with     x(a) = xa
  19. %      y' = f(t,x,y)               y(a) = ya
  20.  
  21. % The formula for f(t,x,y) and is used to form the
  22. %    vector function F(t,Y)    which is stored in  fn.m 
  23.  
  24. % function W = fn(t,Z)
  25. % x = Z(1);  y = Z(2);
  26. % W = [y, -40*sin(x)];
  27.  
  28. delete fn.m
  29. diary fn.m; disp('function W = fn(t,Z)');...
  30.             disp('x = Z(1);  y = Z(2);');...
  31.             disp('W = [y, -40*sin(x)];');...
  32. diary off;
  33.  
  34. fn(0,[0 0]); % Remark.  fn.m and rks4 are used in Algorithm 9.H
  35. pause           % Press any key to continue.
  36.  
  37. clc;
  38. %    Place the endpoints of [a,b] in  a  and  b.
  39.  
  40. %    Place the initial value  Z(a) = [xa ya] in  Za.
  41.  
  42. %    Place the number of subintervals in  m.
  43.  
  44. a  = 0;
  45. b  = 2;
  46. Za = [0.3  0.0];
  47. m  = 40;
  48. [T,Z] = rks4('fn',a,b,Za,m);
  49. P = [T;Z(:,1)']';
  50. points = P(1:4:length(P),:);
  51.  
  52. pause    % Press any key to see the list of solution points.
  53.  
  54. clc;, clg;
  55. X = Z(:,1);
  56. a = min(T)-0.3;
  57. b = max(T)+0.3;
  58. c = min(X)-0.3;
  59. d = max(X)+0.3;
  60. axis([0 2 -0.35 0.35]);...
  61. plot(T,X,'g');...
  62. hold on;...
  63. plot([0 2],[0 0],'r',[0 0],[-0.35 0.35],'r');...
  64. xlabel('t');...
  65. ylabel('x');...
  66. title('Runge-Kutta solution to x`` = f(t,x,x`)');...
  67. grid;...
  68. axis;...
  69. hold off;...
  70. shg; pause    % Press any key to continue.
  71.  
  72. Mx1 = 'Runge-Kutta solution to x`` = f(t,x,x`).';
  73. Mx2 = '     t(k)               x(k)';
  74. clc,echo off,diary output,...
  75. disp(''),disp(Mx1),...
  76. disp(''),disp(Mx2),disp(points),diary off,echo on
  77.